DALS010-Matrix Algebra(矩阵代数)02-Matrix Notation

前言

这一部分是《Data Analysis for the life sciences》的第4章：矩阵代数的第2小节，这一部分的主要内容涉及矩阵的表示方法。

线性模型语言

线性代数的符号化实际上简化了线性模型的数学描述和操作，以及R中的代码。这一部分的练习主要是通过使用矩阵符号来表示线性模型，并使用它来解最小二乘方程。

解方程组

数学家创建线性代数是为了解以下这样的线性方程组：

$$\begin{array}{r}{a+b+c=6} \\ {3 a-2 b+c=2} \\ {2 a+b-c=1}\end{array}$$

它提供了非常有用的机制来解决这些常见问题。我们将会学习如何使用矩阵代数来表示以及求解这些方程组，如下所示：

$$\left(\begin{array}{ccc}{1} & {1} & {1} \\ {3} & {-2} & {1} \\ {2} & {1} & {-1}\end{array}\right)\left(\begin{array}{l}{a} \\ {b} \\ {c}\end{array}\right)=\left(\begin{array}{l}{6} \\ {2} \\ {1}\end{array}\right) \Rightarrow\left(\begin{array}{l}{a} \\ {b} \\ {c}\end{array}\right)=\left(\begin{array}{ccc}{1} & {1} & {1} \\ {3} & {-2} & {1} \\ {2} & {1} & {-1}\end{array}\right)^{-1}\left(\begin{array}{l}{6} \\ {2} \\ {1}\end{array}\right)$$

这一部分就是要解释上面的这些符号是什么意思，以及在统计学中如何用这些符号来计算线性方程。

向量，矩阵与标量

线性代数中的基本元素包括向量(vector)，矩阵(matrix)和标题(scalar)。

向量

在自由落体运行，父子身高以及小鼠体重实验中，与这些数据有关的随机变量我们通常用$Y_{1}, \ldots, Y_{n}$来表示，我们可以将它们视为一个向量(Vector)，其实R中就是这么干的，如下所示：

library(UsingR)
y=father.son$fheight
head(y)

如下所示：

> head(y)
[1] 65.04851 63.25094 64.95532 65.75250 61.13723 63.02254

矩阵

在数学中，我们可以只用一个符号来表示向量，我们通常用一个加粗的黑体字母，即$\mathbf{Y}$来表示向量，从而将它与$Y$区分开来，如下所示：

$$\mathbf{Y}=\left(\begin{array}{c}{Y_{1}} \\ {Y_{2}} \\ {\vdots} \\ {Y_{N}}\end{array}\right)$$

因此，我们知道，一个数据向量的默认就是$N \times 1$维，而不是$1 \times N$维。这里不使用加粗字母是因为，在文中通常会告诉我们这是一个矩阵，而不是一个向量。类似的，我们可以使用数学符合来表示协方差(covariates)或预测因子(predictors)，如果我们有两个预测因子，那么就可以按以下形式表示：

$$\mathrm{X}_{1}=\left(\begin{array}{c}{x_{1,1}} \\ {\vdots} \\ {x_{N, 1}}\end{array}\right) \text { and } \mathrm{X}_{2}=\left(\begin{array}{c}{x_{1,2}} \\ {\vdots} \\ {x_{N, 2}}\end{array}\right)$$

在自由落体案例中，$x_{1,1}=t_{i}$，$x_{i, 1}=t_{i}^{2}$，其中，$t_{i}$表示的是第i次的观测值，这里再强调一次，向量可以被视为$N \times 1$矩阵。因此上面的两个预测因子可以转换为矩阵，如下所示：

$$\mathbf{X}=\left[\mathbf{X}_{1} \mathbf{X}_{2}\right]=\left(\begin{array}{cc}{x_{1,1}} & {x_{1,2}} \\ {\vdots} \\ {x_{N, 1}} & {x_{N, 2}}\end{array}\right)$$

此时，这个矩阵就是$N \times 2$维，在R中，我们可以创建这个矩阵，如下所示：

n <- 25
tt <- seq(0,3.4,len=n) ##time in secs, t is a base function
X <- cbind(X1=tt,X2=tt^2)
head(X)
dim(X)

如下所示：

> head(X)
            X1         X2
[1,] 0.0000000 0.00000000
[2,] 0.1416667 0.02006944
[3,] 0.2833333 0.08027778
[4,] 0.4250000 0.18062500
[5,] 0.5666667 0.32111111
[6,] 0.7083333 0.50173611
> dim(X)
[1] 25  2

我们也可以使用这些符号来表示任意$N \times p$维的矩阵，如下所示：

$$\mathbf{X}=\left(\begin{array}{ccc}{x_{1,1}} & {\dots} & {x_{1, p}} \\ {x_{2,1}} & {\dots} & {x_{2, p}} \\ {} & {\vdots} \\ {x_{N, 1}} & {\dots} & {x_{N, p}}\end{array}\right)$$

在R中，可以使用matrix()函数来创建矩阵，如下所示：

N <- 100; p <- 5
X <- matrix(1:(N*p),N,p)
head(X)
dim(X)

结果如下所示：

> head(X)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1  101  201  301  401
[2,]    2  102  202  302  402
[3,]    3  103  203  303  403
[4,]    4  104  204  304  404
[5,]    5  105  205  305  405
[6,]    6  106  206  306  406
> dim(X)
[1] 100   5

默认情况下，在R中，矩阵是按照列的顺序进行填充的，如果加上参数byrow=TRUE，则矩阵会按照行进行填充，如下所示：

N <- 100; p <- 5
X <- matrix(1:(N*p),N,p,byrow=TRUE)
head(X)

结果如下所示：

> head(X)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    6    7    8    9   10
[3,]   11   12   13   14   15
[4,]   16   17   18   19   20
[5,]   21   22   23   24   25
[6,]   26   27   28   29   30

标量

标量仅仅是一个数字，通常使用小写字母来表示标量，并且不加粗。

练习

P157

数学运算

前面我们提到了一个方程组，如下所示：

$$\begin{aligned} a+b+c &=6 \\ 3 a-2 b+c &=2 \\ 2 a+b-c &=1 \end{aligned}$$

如果用矩阵代数来表示，就是下面的形式：

$$\left(\begin{array}{ccc}{1} & {1} & {1} \\ {3} & {-2} & {1} \\ {2} & {1} & {-1}\end{array}\right)\left(\begin{array}{l}{a} \\ {b} \\ {c}\end{array}\right)=\left(\begin{array}{l}{6} \\ {2} \\ {1}\end{array}\right) \Rightarrow\left(\begin{array}{l}{a} \\ {b} \\ {c}\end{array}\right)=\left(\begin{array}{ccc}{1} & {1} & {1} \\ {3} & {-2} & {1} \\ {2} & {1} & {-1}\end{array}\right)^{-1}\left(\begin{array}{l}{6} \\ {2} \\ {1}\end{array}\right)$$

矩阵与标量相乘

这是矩阵运算中最简单的操作，一个矩阵$\mathbf{X}$与一个标量$a$相乘，矩阵的每一个元素都与这个标量相乘，如下所示：

$$\mathbf{X}=\left(\begin{array}{ccc}{x_{1,1}} & {\dots} & {x_{1, p}} \\ {x_{2,1}} & {\dots} & {x_{2, p}} \\ {} & {\vdots} & {} \\ {x_{N, 1}} & {\dots} & {x_{N, p}}\end{array}\right) \Rightarrow a \mathbf{X}=\left(\begin{array}{ccc}{a x_{1,1}} & {\dots} & {a x_{1, p}} \\ {a x_{2,1}} & {\dots} & {a x_{2, p}} \\ {} & {\vdots} \\ {a_{N, 1}} & {\dots} & {a x_{N, p}}\end{array}\right)$$

R中的运算会自动遵循这个规则，如下所示：

X <- matrix(1:12,4,3)
print(X)
a <- 2
print(a*X)

结果如下所示：

> X <- matrix(1:12,4,3)
> print(X)
     [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12
> a <- 2
> print(a*X)
     [,1] [,2] [,3]
[1,]    2   10   18
[2,]    4   12   20
[3,]    6   14   22
[4,]    8   16   24

转置(Transpose)

矩阵的转置就是将矩阵的行与列颠倒，通常使用$T$这个符号表示矩阵的转置运算，如下所示：

$$\mathbf{X}=\left(\begin{array}{ccc}{x_{1,1}} & {\ldots} & {x_{1, p}} \\ {x_{2,1}} & {\ldots} & {x_{2, p}} \\ {} & {\vdots} & {} \\ {x_{N, 1}} & {\ldots} & {x_{N, p}}\end{array}\right) \Rightarrow \mathbf{X}^{\top}=\left(\begin{array}{ccc}{x_{1,1}} & {\ldots} & {x_{p, 1}} \\ {x_{1,2}} & {\ldots} & {x_{p, 2}} \\ {} & {\vdots} & {} \\ {x_{1, N}} & {\ldots} & {x_{p, N}}\end{array}\right)$$

在R中，直接使用t()函数就行，如下所示：

X <- matrix(1:12,4,3)
X
t(X)

如下所示：

> X <- matrix(1:12,4,3)
> X
     [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12
> t(X)
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    5    6    7    8
[3,]    9   10   11   12

矩阵的相乘

还以开始三元一次方程组为例说明一下，如下所示：

$$\begin{array}{r}{a+b+c=6} \\ {3 a-2 b+c=2} \\ {2 a+b-c=1}\end{array}$$

两个矩阵相乘的运算法则是，第1个矩阵的行乘以第2个矩阵的列（第1个矩阵的列数与第2个矩阵的行数相同）。由于第2个矩阵只有1列，因此相乘的结果如下所示：

$$\left(\begin{array}{ccc}{1} & {1} & {1} \\ {3} & {-2} & {1} \\ {2} & {1} & {-1}\end{array}\right)\left(\begin{array}{l}{a} \\ {b} \\ {c}\end{array}\right)=\left(\begin{array}{c}{a+b+c} \\ {3 a-2 b+c} \\ {2 a+b-c}\end{array}\right)$$

下面我们在R中计算一下，我们来看一下abc=c(3,2,1)是否是上述方程组的一个解，计算结果如下所示：

X <- matrix(c(1,3,2,1,-2,1,1,1,-1),3,3)
abc <- c(3,2,1) #use as an example
rbind( sum(X[1,]*abc), sum(X[2,]*abc), sum(X[3,]%*%abc))

计算结果如下所示：

> rbind( sum(X[1,]*abc), sum(X[2,]*abc), sum(X[3,]%*%abc))
     [,1]
[1,]    6
[2,]    6
[3,]    7

使用%*%可以进行矩阵相乘，这种方法很简单，如下所示：

X%*%abc

结果如下所示：

> X%*%abc
     [,1]
[1,]    6
[2,]    6
[3,]    7

从上面结果我们可以知道，c(3,2,1)不是方程组的解。为了求解方程组，我们需要逆转(inverse)左边的矩阵（后面要学习的内容）。在这里，我们定义矩阵$A$与$X$相乘的一般规则，如下所示：

$$\mathbf{A X}=\left(\begin{array}{cccc}{a_{1,1}} & {a_{1,2}} & {\dots} & {a_{1, N}} \\ {a_{2,1}} & {a_{2,2}} & {\dots} & {a_{2, N}} \\ {} & {} & {\vdots} \\ {a_{M, 1}} & {a_{M, 2}} & {\dots} & {a_{M, N}}\end{array}\right)\left(\begin{array}{ccc}{x_{1,1}} & {\dots} & {x_{1, p}} \\ {x_{2,1}} & {\dots} & {x_{2, p}} \\ {} & {\vdots} \\ {x_{N, 1}} & {\dots} & {x_{N, p}}\end{array}\right)\\=\left(\begin{array}{ccc}{\sum_{i=1}^{N} a_{1, i} x_{i, 1}} & {\ldots} & {\sum_{i=1}^{N} a_{1, i} x_{i, p}} \\ {} & {\vdots} & {} \\ {\sum_{i=1}^{N} a_{M, i} x_{i, 1}} & {\ldots} & {\sum_{i=1}^{N} a_{M, i} x_{i, p}}\end{array}\right)$$

只有当矩阵$A$的列数与矩阵$X$的行数相等时，才能得到两个矩阵相乘的结果。

单位矩阵(The identity matrix)

在矩阵代数中，单位矩阵的意义与1类似，一个矩阵乘以单位矩阵，还是其本身，单位矩阵的定义如下所示：

$$\mathbf{I}=\left(\begin{array}{cccccc}{1} & {0} & {0} & {\ldots} & {0} & {0} \\ {0} & {1} & {0} & {\ldots} & {0} & {0} \\ {0} & {0} & {1} & {\ldots} & {0} & {0} \\ {\vdots} & {\vdots} & {\vdots} & {\ddots} & {\vdots} & {\vdots} \\ {0} & {0} & {0} & {\ldots} & {1} & {0} \\ {0} & {0} & {0} & {\ldots} & {0} & {1}\end{array}\right)$$

单位矩阵的行与列数目相等，从左上解到右下解的元素都为1，其余都为0。现在我们看一下一个矩阵与单位矩阵的相乘结果，如下所示：

$$\mathbf{X I}=\left(\begin{array}{ccc}{x_{1,1}} & {\dots} & {x_{1, p}} \\ {} & {\vdots} & {} \\ {x_{N, 1}} & {\dots} & {x_{N, p}}\end{array}\right)\left(\begin{array}{cccccc}{1} & {0} & {0} & {\dots} & {0} & {0} \\ {0} & {1} & {0} & {\dots} & {0} & {0} \\ {0} & {0} & {1} & {\dots} & {0} & {0} \\ {} & {} & {} & {\vdots} & {} \\ {0} & {0} & {0} & {\dots} & {1} & {0} \\ {0} & {0} & {0} & {\dots} & {0} & {1}\end{array}\right)=\left(\begin{array}{ccc}{x_{1,1}} & {\dots} & {x_{1, p}} \\ {} & {\vdots} & {} \\ {x_{N, 1}} & {\dots} & {x_{N, p}}\end{array}\right)$$

R中可以通过diag()函数来生成单位矩阵，如下所示：

n <- 5 #pick dimensions
diag(n)

结果如下所示：

> n <- 5 #pick dimensions
> diag(n)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    0    1    0    0    0
[3,]    0    0    1    0    0
[4,]    0    0    0    1    0
[5,]    0    0    0    0    1

矩阵的逆转(inverse)

矩阵$X$的逆转可以表示为$X^{-1}$，逆转后的矩阵与原矩阵相乘，会得到单位矩阵，即$X^{-1} X=I$，但是，并非所有的矩阵都可以逆转（逆转后的矩阵称为逆矩阵）。

当我们计算一个线性模型时，就要用到逆矩阵。在R中可以使用solve()函数来进行运算，solve()函数就是计算一个矩阵的逆矩阵，如下所示：

X <- matrix(c(1,3,2,1,-2,1,1,1,-1),3,3)
y <- matrix(c(6,2,1),3,1)
solve(X)
solve(X)%*%y #equivalent to solve(X,y)

结果如下所示：

> solve(X)
           [,1]        [,2]       [,3]
[1,] 0.07692308  0.15384615  0.2307692
[2,] 0.38461538 -0.23076923  0.1538462
[3,] 0.53846154  0.07692308 -0.3846154
> solve(X)%*%y #equivalent to solve(X,y)
     [,1]
[1,]    1
[2,]    2
[3,]    3

练习

P163

矩阵案例

前面内容是关于矩阵的简单案例，以及它们在数据分析过程中的重要作用。我们最终要达到的目标是：如何使用矩阵代数符号来描述线性模型以及求解最小二乘问题。

均值

为了计算样本的均值和方差，例如我们使用$\overline{Y}=\frac{1}{N} Y_{i}$来计算均值，使用$\operatorname{var}(Y)=\frac{1}{N} \sum_{i=1}^{N}\left(Y_{i}-\overline{Y}\right)^{2}$来计算方差。我们也可以使用矩阵乘法来计睡这些结果。

第一步：定义一个$N\times1$矩阵，它只有1列，如下所示：

$$A=\left(\begin{array}{l}{1} \\ {1} \\ {\vdots} \\ {1}\end{array}\right)$$

我们可以使用下面的公式来计算均值：

$$\frac{1}{N} \mathbf{A}^{\top} Y=\frac{1}{N}\left(\begin{array}{llll}{1} & {1} & {\ldots} & {1}\end{array}\right)\left(\begin{array}{l}{Y_{1}} \\ {Y_{2}} \\ {\vdots} \\ {Y_{N}}\end{array}\right)=\frac{1}{N} \sum_{i=1}^{N} Y_{i}=\overline{Y}$$

从上面的公式我们可以看到，里面乘了一个标量，即$1/N$，在R中，可以使用%*%来实现，如下所示：

library(UsingR)
y <- father.son$sheight
print(mean(y))
N <- length(y)
Y<- matrix(y,N,1)
A <- matrix(1,N,1)
barY=t(A)%*%Y / N
print(barY)

结果如下所示：

> library(UsingR)
> y <- father.son$sheight
> print(mean(y))
[1] 68.68407
> N <- length(y)
> Y<- matrix(y,N,1)
> A <- matrix(1,N,1)
> barY=t(A)%*%Y / N
> print(barY)
         [,1]
[1,] 68.68407

方差

在统计学中，常常会用到矩阵转置后的相乘，而在R中，可以直接进行运算，如下所示：

barY=crossprod(A,Y) / N
print(barY)
barY2 <- t(A)%*%Y/N
print(barY2)

计算结果如下所示：

> print(barY)
         [,1]
[1,] 68.68407
> barY2 <- t(A)%*%Y/N
> print(barY2)
         [,1]
[1,] 68.68407

在计算方差时，公式如下所示：

$$\mathbf{r} \equiv\left(\begin{array}{c}{Y_{1}-\overline{Y}} \\ {\vdots} \\ {Y_{N}-\overline{Y}}\end{array}\right), \frac{1}{N} \mathbf{r}^{\top} \mathbf{r}=\frac{1}{N} \sum_{i=1}^{N}\left(Y_{i}-\overline{Y}\right)^{2}$$

在R中，如果crossprod()函数只有一个矩阵参数，那么计算的就是$r^{\top} r$，可以简单地写为以下形式：

r <- y - barY
crossprod(r)/N
t(r)%*%r/N

结果如下所示：

> crossprod(r)/N
         [,1]
[1,] 7.915196
> t(r)%*%r/N
         [,1]
[1,] 7.915196

上述计算过程也等于以下形式：

library(rafalib)
popvar(y)

结果如下所示：

> popvar(y)
[1] 7.915196

线性模型

现在我们应用一下线性模型，还以Galton的案例为例说明一下，我们先定义以下矩阵：

$$\mathbf{Y}=\left(\begin{array}{c}{Y_{1}} \\ {Y_{2}} \\ {\vdots} \\ {Y_{N}}\end{array}\right), \mathbf{X}=\left(\begin{array}{cc}{1} & {x_{1}} \\ {1} & {x_{2}} \\ {\vdots} & {} \\ {1} & {x_{N}}\end{array}\right), \beta=\left(\begin{array}{c}{\beta_{0}} \\ {\beta_{1}}\end{array}\right) \text { and } \varepsilon=\left(\begin{array}{c}{\varepsilon_{1}} \\ {\varepsilon_{2}} \\ {\vdots} \\ {\varepsilon_{N}}\end{array}\right)$$

然后写出线性模型，如下所示：

$$Y_{i}=\beta_{0}+\beta_{1} x_{i}+\varepsilon, i=1, \ldots, N$$

它等于以下形式：

$$\left(\begin{array}{c}{Y_{1}} \\ {Y_{2}} \\ {\vdots} \\ {Y_{N}}\end{array}\right)=\left(\begin{array}{cc}{1} & {x_{1}} \\ {1} & {x_{2}} \\ {\vdots} & {} \\ {1} & {x_{N}}\end{array}\right)\left(\begin{array}{c}{\beta_{0}} \\ {\beta_{1}}\end{array}\right)+\left(\begin{array}{c}{\varepsilon_{1}} \\ {\varepsilon_{2}} \\ {\vdots} \\ {\varepsilon_{N}}\end{array}\right)$$

简化一下，就是下面的形式：

$$\mathrm{Y}=\mathrm{X} \beta+\varepsilon$$

那么最小二乘方程就变得比较简单，如下所示：

$$(\mathrm{Y}-\mathrm{X} \boldsymbol{\beta})^{\top}(\mathrm{Y}-\mathrm{X} \boldsymbol{\beta})$$

现在我们计算的目标就是，找到$\beta$的值，使得上述方程的值最小，我们可以使用微积分的手段来进行计算。

线性模型的微积分计算

在使用矩阵符号来计算偏微分方程(partial derivative equation)时，需要遵循几个规则。通过当微分方程为0时，计算出$\beta$，这个就是我们要求的解。在这里，我们上述的微分方程是以下形式：

$$\begin{array}{c}{2 \mathbf{X}^{\top}(\mathbf{Y}-\mathbf{X} \hat{\boldsymbol{\beta}})=0} \\ {\mathbf{X}^{\top} \mathbf{X} \hat{\boldsymbol{\beta}}=\mathbf{X}^{\top} \mathbf{Y}} \\ {\hat{\boldsymbol{\beta}}=\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \mathbf{Y}}\end{array}$$

对于求出的解，也就是$\beta$，通常在上面加一个帽子结构，即$\hat{\beta}$，这个表示是对真实的$\beta$值的估计。这里我们需要记住的是，最小二乘类似于一个幂运算，这个公式类似于$f(x)^{2}$的导数2$f(x) f^{\prime}(x)$。

在R中计算LSE

代码如下所示：

library(UsingR)
x=father.son$fheight
y=father.son$sheight
X <- cbind(1,x)
betahat <- solve( t(X) %*% X ) %*% t(X) %*% y
###or
betahat <- solve( crossprod(X) ) %*% crossprod( X, y )

结果如下所示：

> betahat
       [,1]
  33.886604
x  0.514093

通过计算估计的$\hat{\beta_{0}}+\hat{\beta_{1}}x$就能计算出相应的$x$的值，如下所示：

newx <- seq(min(x),max(x),len=100)
X <- cbind(1,newx)
fitted <- X%*%betahat
plot(x,y,xlab="Father's height",ylab="Son's height")
lines(newx,fitted,col=2)

其中，$\hat{\boldsymbol{\beta}}=\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \mathbf{Y}$是数据分析中最常用的公式之一，它的一大优势就是可以应用于多种情况，例如在自由落体运算中：

set.seed(1)
g <- 9.8 #meters per second
n <- 25
tt <- seq(0,3.4,len=n) #time in secs, t is a base function
d <- 56.67 - 0.5*g*tt^2 + rnorm(n,sd=1)

我们几乎用了与前面相同的代码，如下所示：

X <- cbind(1,tt,tt^2)
y <- d
betahat <- solve(crossprod(X))%*%crossprod(X,y)
newtt <- seq(min(tt),max(tt),len=100)
X <- cbind(1,newtt,newtt^2)
fitted <- X%*%betahat
plot(tt,y,xlab="Time",ylab="Height")
lines(newtt,fitted,col=2)

最终的估计值为：

> betahat
         [,1]
   56.5317368
tt  0.5013565
   -5.0386455

lm()函数

现在我们使用lm()函数来计算一下自由落体运动，如下所示：

X <- cbind(tt,tt^2)
fit=lm(y~X)
summary(fit)

结果如下所示：

> summary(fit)

Call:
lm(formula = y ~ X)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.5295 -0.4882  0.2537  0.6560  1.5455 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  56.5317     0.5451 103.701   <2e-16 ***
Xtt           0.5014     0.7426   0.675    0.507    
X            -5.0386     0.2110 -23.884   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9822 on 22 degrees of freedom
Multiple R-squared:  0.9973,	Adjusted R-squared:  0.997 
F-statistic:  4025 on 2 and 22 DF,  p-value: < 2.2e-16

这个结果与前面的计算结果一致。

总结

在这一部分中，我们学习了如何使用线性代数来描述线性模型。随后我们研究了几个案例，其中的一些是与实验设计有有。我们还演示了最小二乘估计。但是，我们需要记住，由于y是一个随机变量，这些估计也是随机的。在后面的章节中，我们将学习如何计算这些随机变量，以及随机变量的标准误以及统计推断。

练习

P171